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We  study  the  dynamical  behavior  of  one-dimensional  arrays  of  ~100  pm  diameter  aqueous  droplets 
containing  the  oscillatory  Belousov-Zhabotinsky  (BZ)  reaction,  separated  by  narrow  gaps  of  a  fluorinated 
oil.  In  this  closed  system,  the  malonic  acid  concentration  decreases  as  the  reaction  proceeds.  Starting 
with  a  low  initial  malonic  acid  concentration,  we  observe  a  series  of  attractors  as  a  function  of  time  in  the 
following  order:  anti-phase  attractors;  in-phase  attractors,  which  evolve  into  traveling  waves;  and  mixed 
modes  that  contain  either  regions  of  in-phase  droplets  separated  by  anti-phase  oscillators,  or  in-phase 
oscillators  combined  with  non-oscillatory  droplets.  Most  of  the  observations  are  consistent  with  numerical 
chemical  models  of  the  BZ  reaction  in  which  components  that  participate  in  the  excitatory  (bromine 
dioxide  and  bromous  acid)  and  inhibitory  (bromine)  pathways  diffuse  between  the  droplets.  The  models 
are  used  to  quantitatively  assess  the  inter-drop  coupling  strength  as  a  function  of  drop  separation,  drop 
size  and  malonic  acid  concentration.  To  experimentally  establish  the  mechanism  of  excitatory  coupling 
between  the  BZ  droplets,  we  verify  the  transport  through  the  fluorinated  oil  of  chlorine  dioxide  and 
several  weak  acids,  including  malonic  acid. 


1.  Introduction 

Recent  studies1-8  of  chemically  coupled  microdroplets  containing 
the  Belousov-Zhabotinsky  (BZ)  oscillating  chemical  reaction 
demonstrate  that  such  emulsion  systems  are  interesting  experi¬ 
mental  models  for  understanding  coupled  nonlinear  chemical 
oscillators.  These  systems  show  collective  behavior  belonging  to 
the  same  class  of  phenomena  exhibited  by  numerous  biological 
systems,  ranging  from  single  cell  organisms  to  neuronal  circuits. 
Specifically,  the  system  studied  here  consists  of  drops  of  BZ 
solution  of  ~  100  pm  diameter  immersed  in  oil.  The  approxi¬ 
mately  one  nanoliter  drops  are  stabilized  by  a  surfactant  and 
enclosed  in  a  cylindrical  capillary  of  inner  diameter  of  100  pm, 
forming  a  linear  one-dimensional  (l-D)  array  of  about  one 
hundred  BZ  drops,  as  illustrated  in  Fig.  1.  The  BZ  reactants, 
intermediates  and  products  partition  into  the  oil  phase  to 
different  degrees  and  can  diffuse  between  the  drops,  causing 
the  coupling.  In  contrast  to  previous  work,1’3,5,8  in  which  we 
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studied  drops  in  1-D  and  2-D  geometries  and  observed  a  rich 
variety  of  dynamical  attractor  states  ranging  from  all  the  drops 
oscillating  with  fixed  phase  shifts  between  neighboring  drops, 
to  stationary  Turing  states,  as  well  as  mixed  states  in  which 
some  drops  were  stationary  and  others  oscillated,  in  this  paper  we 
study  1-D  drop  dynamics  at  extremes  of  high  and  low  malonic 
acid  concentration  ([MA]).  A  previous  paper2  is  devoted  to  theore¬ 
tical  modeling  of  the  observations  presented  here.  In  the  limit  of 
low  [MA],  we  also  observe  in-phase  oscillations  and  traveling 
waves.  All  our  models2  suggest  that  the  in-phase  behavior  at 
low  [MA]  is  due  to  inter-drop  coupling  via  excitatory  species. 
Additionally,  we  observe  mixed  mode  patterns  in  which  some 
drops  oscillate  and  others  do  not,  as  well  as  “local  in-phase” 
behavior,  in  which  pairs  or  small  clusters  of  neighboring  droplets 
are  in  phase  with  each  other,  but  are  out-of-phase  with  neighboring 
clusters.  This  behavior  is  consistent  with  multistable  attractors 
resulting  from  a  combination  of  inhibitory  and  excitatory 
coupling.2  In  the  high  [MA]  limit,  we  also  observe  in-phase 
oscillations  and  traveling  waves.  However,  our  models2  do  not 
predict  in-phase  behavior  at  high  [MA],  indicating  that  our 
chemical  mechanism  of  inter-drop  coupling  is  incomplete. 

In  all  cases  we  compare  experiments  on  closed  systems  with 
theory  for  open  systems  in  which  the  concentrations  of  malonic 
acid,  sulfuric  acid  and  bromate  ion  are  assumed  to  remain  constant. 
This  enables  us  to  identify  stable  attractors  and  to  avoid  modeling 
the  poorly  understood  details  of  malonic  acid  consumption. 
As  long  as  the  fractional  change  in  the  concentrations  of  these 
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Fig.  1  BZ  oscillations  in  a  1-D  closed  array  of  drops,  (a)  Photograph  of 
section  of  a  hydrophobized  cylindrical  capillary  of  100  pm  inner  diameter 
filled  with  drops  containing  the  Belousov-Zhabotinsky  reaction.  The  drops 
are  separated  by  fluorinated  oil.  (b)  Magnified  section  of  an  unhydropho- 
bized  capillary  containing  BZ  drops  (red  arrows)  separated  by  oil  (green 
arrows).  Droplet  diameter,  a  =  120  |im.  Oil  gap,  b  =  40  |im.  (c)  Space-time 
plot  of  seven  drops  shown  in  Fig.  lb.  Three  arrows  indicate  the  oil  gaps 
between  drops,  which  appear  as  bright  lines  running  parallel  to  the  time 
axis,  because  the  light  transmitted  through  the  oil  gap  is  constant  in  time. 
Bright  lines  parallel  to  the  space  axis  correspond  to  the  oxidized  state  of 
the  BZ  reaction,  (d)  Period  vs.  time  for  the  seven  drops  shown  in  the  space- 
time  plot  of  Fig.  lc.  Initially,  [MA]  =  40  mM.  (e)  Transmitted  light  intensity 
through  drop  6  (from  left)  from  Fig.  10a  as  a  function  of  time.  The  intensity 
profile  as  a  function  of  time  for  the  2nd,  15th,  and  22nd  oscillations, 
(f)  Calculated  ferroin  concentration  vs.  time  for  different  concentrations 
of  malonic  acid,  with  [Br03“]  =  600  mM  and  [H2S04]  =  80  mM,  a  =  b  = 
150  |!m.  (g)  Calculated  oscillation  period  vs.  malonic  acid  concentration  for 
a  single  BZ  droplet,  using  default  chemical  conditions:  [Br03_]  =  300  mM 
and  [H2S04]  =  80  mM. 


three  chemicals  is  small,  the  experiments  on  the  closed  systems 
and  our  simulations  of  the  open  systems  agree  well.  Here  we 
consider  systems  in  which  the  concentration  of  MA  is  smaller 
than  the  sulfuric  acid  and  bromate  ion  concentrations,  so  the 
consumption  of  MA  represents  the  greatest  difference  between 
the  open  and  closed  systems.  To  compare  theory  and  experiment, 
we  prepare  a  series  of  BZ  emulsions  as  a  function  of  [MA]  and 
compare  the  initial  behavior  of  each  of  these  emulsions,  which  are 
closed  systems,  with  the  corresponding  theoretical  open  systems. 
Pursuing  the  same  logic,  we  reason  that  the  time  dependence  of  a 
single  closed  system  can  be  qualitatively  mapped  onto  a  series  of 
open  systems  with  decreasing  [MA]. 

Theoretical  models  of  coupled  BZ  oscillators  are  necessarily 
simplified,  not  only  because  of  the  complexity  of  the  chemical 
kinetics  of  the  reaction,9  but  also  because  the  degree  to  which 
the  chemicals  partition  into  the  oil  and  their  reactivity  within  the 
oil  are  poorly  known.  However,  because  the  kinetics  within  the 
microdroplets  are  the  same  as  in  the  macroscopic  BZ  reaction, 
one  can  make  use  of  well-established  detailed  models10  as  a 
starting  point.3,5 

We  model  the  coupled  BZ  drops  at  three  levels  of  approxi¬ 
mation.  The  most  realistic  consists  of  a  finite  element  model 
(COMSOL®)  of  the  reaction-diffusion  equation  in  1-D,  2-D  and 
3-D,  where  each  drop  is  modeled  as  a  point,  disk  or  sphere, 
respectively.  Only  minor  differences  are  observed  between  1-D, 
2-D  and  3-D  drops  arranged  periodically  in  a  straight  line.  The 
capillary  walls  are  treated  as  a  no-flux  boundary  condition.  The 
reaction  term  is  modeled  using  the  7-variable  FKN  model  as 
detailed  in  previous  work.1  Chemicals  are  allowed  to  diffuse 
within  the  drop  and  the  oil  and  to  exchange  between  the  two 
phases  governed  by  a  partition  coefficient,  P,  defined  as  the 
ratio  of  equilibrated  concentrations  of  a  particular  species  in  the 
oil  and  water  phases  (P  =  c0n/cwater).  These  boundary  conditions  at 
the  oil/water  interface  are  enforced  in  COMSOL  using  the  “stiff” 
method11  to  have  to  right  partition  coefficients.  Each  chemical 
species  has  the  same  diffusion  constant  in  the  oil  and  the  water. 
No  chemical  reactions  are  allowed  to  occur  in  the  oil,  which  we 
find  experimentally  is  not  the  case.  However,  we  argue  that  this 
simplification  is  justified  for  most  of  the  cases  studied  here. 
Models  of  intermediate  complexity  approximate  each  spatially 
extended  drop  as  a  7-variable  FKN  oscillator  confined  to  a  single 
point.  The  BZ  drops  are  diffusively  coupled  to  their  nearest 
neighbors.  This  class  of  models,  consisting  of  a  ring  of  diffusively 
coupled  point-like  chemical  reactors,  was  first  introduced 
by  Turing.12  We  numerically  solve  the  model  using  Matlab® 
as  described  in  an  earlier  computational  study.1  The  diffusion 
across  a  100  pm  sized  droplet  takes  very  little  time  compared 
to  the  oscillation  period,  so  the  droplets  are  oscillating  homo¬ 
geneously.  Thus  the  point  model  is  justified.  The  point  model 
and  finite-element  model  give  similar  results  for  1-D  systems. 
Finally,  at  the  most  abstract  level,  we  construct  a  single  variable 
phase  oscillator  model,13  in  which  the  phase  coupling  function 
is  derived  from  the  point  oscillator  model. 

In  our  numerical  investigations  of  open  systems,  steady  state 
dynamical  patterns  are  identified  as  attractors.  Technically,  no 
attractors  exist  in  the  closed  system,  because  the  only  steady 
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state  is  equilibrium.  However,  as  we  previously  demonstrated, 
our  closed  system  gives  rise  to  “pseudo-attractors”  that  main¬ 
tain  an  oscillatory  pattern  for  many  periods  and  have  the  same 
dynamic  behavior  as  the  true  attractors  found  in  simulations  of 
the  corresponding  open  system.1 

We  also  experimentally  explore  the  source  of  stability  of  the 
in-phase  attractors  in  our  system;  experiments  to  measure  the 
transport  of  several  BZ-related  species  through  our  fluorinated 
oil  suggest  that  excitatory  coupling  takes  place  via  exchange  of 
bromine  dioxide  radical  (Br02*)  and  bromous  acid  (HBr02) 
among  droplets,  while  the  inhibitory  coupling  is  conveyed  by 
bromine  (Br2). 

2.  Experimental 

Arrays  of  water  droplets  in  oil  are  generated  using  flow  focusing 
microfluidic  technology  as  described  in  our  previous  work:1 
A  microfluidic  chip  fabricated  from  polydimethylsiloxane  (PDMS) 
allows  the  injection  of  equal  amounts  of  different  aqueous  streams 
into  a  central  mixing  channel.  This  channel  ends  in  a  nozzle  of 
50  pm  diameter,  where  streams  of  fluorinated  oil  coming  perpen¬ 
dicularly  from  both  sides  produce  droplets  by  flow  focusing.  After 
the  nozzle,  the  channel  widens  to  100  pm  diameter  with  a  square 
cross  section.  Drops  of  this  volume  oscillate  at  most  a  few  times  if 
stored  in  a  PDMS  channel.  Previously  we  argued  this  is  due  to 
permeation  of  bromine  into  the  PDMS.  If  the  ratio  of  the  volume  of 
the  PDMS  chip  to  the  volume  of  the  BZ  solution  is  large,  then  the 
oscillating  reaction  ceases.3  Therefore  we  export  the  drops  from 
the  PDMS  chip  to  a  hydrophobized  glass  capillary.  In  standard  soft 
lithography  microfluidics,14  three  of  the  four  walls  of  the  micro¬ 
fluidic  channel  are  made  from  PDMS,  which  is  then  sealed  to  a 
glass  slide  to  form  the  fourth  wall.  Instead,  we  use  a  5  mm  flat  slab 
of  PDMS  to  seal  the  channel  and  cut  the  assembled  device 
perpendicular  to  the  channel.  Having  four  soft  walls  allows  us  to 
insert  a  cylindrical  hydrophobized  capillary  of  100  pm  inner 
diameter  (170  pm  outer  diameter)  directly  into  the  100  pm  channel 
emerging  from  the  cut  side  of  the  microfluidic  chip.  Unless 
otherwise  specified,  all  experiments  described  here  are  conducted 
in  capillaries  of  this  size. 

Previous  to  experiments,  capillaries  are  hydrophobized  using  a 
vacuum  chamber  in  which  capillaries  and  a  small  amount  of  liquid 
(tridecafluoro-l,l,2,2-tetrahydrooctyl)trichlorosilane  are  placed.  We 
reduce  the  pressure  in  order  to  evaporate  the  silane,  which  enters 
the  capillaries  via  gaseous  diffusion.  The  trichlorosilane  group  of 
this  molecule  reacts  with  the  oxygen  groups  of  the  silica  on  the 
glass  surface,  covering  it  with  fluorocarbon  tails.15  After  two  hours 
in  the  vacuum  chamber,  which  is  sufficient  for  hydrophobizing 
the  capillaries,  we  remove  the  capillaries.  Air  stops  the  reaction, 
because  oxygen  reacts  with  the  chlorosilane  groups.  The  oil 
used  in  this  work  is  a  perfluorinated  oil,  HFE  7500:  3-ethoxy- 
l,l,l,2,3,4,4,5,5,6,6,6-dodecafluoro-2-trifluoromethyl-hexane  (3M 
Corp.,  St.  Paul,  MN,  USA,  dielectric  constant  =  5.8).  Since  this  is 
the  only  oil  used,  we  refer  it  as  “the  oil”.  An  amphiphilic  PEG- 
PFPE  fluorinated  block  copolymer  dissolved  at  a  concentration  of 
2%  by  volume  in  the  oil  serves  as  a  surfactant  for  aiding  droplet 


formation  and  preventing  droplet  coalescence.  The  surfactant 
was  provided  by  RainDance  Technologies,  Lexington,  MA,  USA 
and  is  similar  to  a  surfactant16  described  in  the  literature.  A 
commercial  version  of  this  surfactant  is  now  available  (RAN 
Biotechnologies). 

To  form  BZ  droplets  in  oil  using  this  chip,  we  inject  malonic 
acid  (Sigma-Aldrich)  and  sulfuric  acid  (Fisher)  at  twice  the  final 
desired  concentrations  in  the  droplet  in  one  aqueous  stream. 
The  second  aqueous  stream  contains  sodium  bromate  (Fisher), 
ferroin  indicator  (Aqua  Solutions)  and  Ru(bipy)32+  (Sigma-Aldrich) 
also  at  twice  the  final  desired  concentrations.  Unless  otherwise 
noted,  [Br03_]  =  300  mM,  [H2S04]  =  80  mM,  [ferroin]  =  3  mM  and 
[Ru(bipy)32+]  =  0.4  mM  in  the  droplets  in  all  our  experiments.  We 
use  ferroin  as  a  BZ  catalyst  and  Ru(bipy)32+  as  a  co-catalyst  and 
photosensitizer.  The  ferroin  contained  in  the  BZ  drops  also  serves 
as  a  dye  indicator,  changing  from  blue  (oxidized  ferroin)  to  red 
(reduced  ferriin)  depending  on  its  oxidation  state.  Capillaries  with 
identical  BZ  droplets  equidistantly  spaced  by  identical  volumes  of 
oil  are  sealed  using  5  minute  epo>^,  attached  to  a  microscope 
slide  and  observed  with  a  CCD  camera  through  a  microscope.  If 
the  capillary  is  inadequately  sealed,  gas  bubbles  evolve,  causing 
the  drops  to  move,  which  disturbs  the  chemical  dynamics  and 
disqualifies  the  experiment. 

The  capillary  is  illuminated  in  transmission  with  a  510  ± 
10  nm  notch  interference  filter  inserted  into  the  light  path  to 
maximize  the  transmission  of  light  when  the  BZ  drops  are  in 
the  oxidized  state.  Visually,  the  BZ  drops  appear  bright  in  the 
ferroin  oxidized  state  (ferriin)  and  dark  in  the  reduced  state. 
The  light  intensity  transmitted  through  the  transparent  oil 
separating  the  drops  is  constantly  bright  during  the  entire 
experiment  because  ferroin  does  not  diffuse  out  of  the  drops 
into  the  oil.  Our  field  of  vision,  in  which  up  to  20  droplets  can 
be  observed,  is  located  at  the  center  of  the  sealed  capillary. 
There  are  at  least  15  additional  droplets  to  each  side  of  this 
central  region.  In  most  cases  the  drops  have  a  few  percent 
variation  in  size  and  spacing.  However,  we  did  not  quantita¬ 
tively  investigate  the  effects  of  variable  drop  size  and  spacing  on 
the  dynamics,  which  will  be  the  subject  of  future  work.  We 
construct  space-time  plots  by  acquiring  the  space  data  at  a  fixed 
time  from  a  single  image  by  extracting  the  intensity  of  trans¬ 
mitted  light  from  a  single  row  of  pixels  along  the  capillary  axis. 
The  time  data  is  acquired  by  repeating  this  process  from  a  set  of 
digital  images  taken  at  constant  (3-5  s)  intervals  over  the  course  of 
an  experiment.  In  some  experiments,  the  droplets  are  initially 
synchronized  in-phase  using  the  photosensitivity  of  the  Ru(bipy)32+ 
co-catalyst  to  suppress  oscillations  by  illuminating  the  capillary 
with  intense  450  nm  light  for  about  10  min.5  The  diameter  of 
the  droplets  is  defined  as  their  width  along  the  capillary  in  the 
axial  direction.  The  distance  between  droplets  is  taken  as  the 
length  of  the  gap  between  the  boundaries  of  adjacent  droplets 
measured  along  the  axis  of  the  capillary. 

Droplets  of  chlorine  dioxide,  C102#,  in  water  are  generated 
using  the  PDMS  microfluidic  chip.  In  one  aqueous  stream,  we 
introduce  sodium  chlorite  at  0.5  M  (NaCl02,  Acros  Organics, 
80%  pure),  pH  =  6,  and  in  the  other  aqueous  stream  we  place 
nitric  acid  in  excess:  4.0  M  (Fisher).  The  acidification  of  NaCl02 
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produces  HC102,  which  undergoes  a  disproportionation  reac¬ 
tion  and  forms  C102*.17  For  the  conditions  used  in  the  droplets, 
spectrophotometric  measurements  in  bulk  reveal  a  maximum 
concentration  of  C102#  of  about  one  tenth  of  the  initial  chlorite 
concentration  approximately  20  min  after  the  mixing  of  the 
reactants.  From  this  observation,  we  can  estimate  the  maximum 
concentration  of  C102*  in  the  droplets  as  ~50  mM.  To  observe 
the  concentration  of  C102#  in  the  microfluidic  drops,  we  place  a 
400  nm  interference  notch  filter  in  the  microscope  illumination 
path.  Although  both  C102#  and  HC102  absorb  at  this  wavelength, 
the  molar  extinction  coefficient  of  C102#  is  at  least  one  order  of 
magnitude  greater  than  that  of  HC102.18 

For  experiments  on  the  transport  of  weak  acids,  thymol  blue, 
iodoacetic  acid  and  chloroacetic  acid  were  obtained  from 
Sigma-Aldrich,  Acros  Organics  and  Fisher,  respectively. 

3.  Results  and  discussion 

This  work  consists  of  four  sections,  (l)  The  quantitative  relation¬ 
ship  between  our  experiments  on  a  closed  BZ  system  and  the 
theory  of  an  open  system.  (2)  The  origin  of  inhibitory  coupling 
between  neighboring  BZ  droplets  and  the  measurement  of  the 
coupling  strength.  (3)  A  phase  model  for  weakly  coupled  BZ 
oscillators.  (4)  The  chemical  mechanisms  of  excitatory  coupling, 
which  gives  various  in-phase  behaviors.  These  points  will  be 
discussed  in  the  following  subsections. 

3.1  Consumption  of  malonic  acid  in  closed  system 

In  Fig.  1  we  present  an  experiment  that  is  representative  of  many 
of  the  phenomena  observed  for  BZ  drops  at  low  malonic  acid. 
Consider  the  complex  behavior  illustrated  in  the  space-time  plot 
of  Fig.  lc.  In  the  beginning,  drops  are  synchronized  to  begin 
oscillating  in-phase,  but  with  time  drift  out-of-phase.  At  the  end, 
some  drops  (e.g,  #4,  #5  and  #6)  switch  back  to  in-phase,  while 
some  others  (e.g.,  #2)  stop  oscillating.  The  oscillation  period 
is  longer  for  drops  that  are  out-of-phase  with  their  neighbors 
than  for  in-phase  drops.  The  period  increases  with  time,  and  the 
oscillation  waveform  of  the  oxidation/reduction  of  ferroin 
changes  with  time.  All  of  these  phenomena  will  be  explained 
as  a  consequence  of  the  consumption  of  malonic  acid,  leading  to 
increased  inter-drop  coupling. 

Proceeding  in  detail,  Fig.  la  is  a  photograph  of  a  100  pm 
inner  diameter  hydrophobized  capillary  filled  with  BZ  drops 
that  are  in  contact  with  each  other.  The  oil  wets  the  capillary, 
and  surface  tension  causes  the  drops  to  adopt  a  cylindrical 
shape  with  hemispheres  on  the  end.  Fig.  lb  shows  120  pm  long 
BZ  drops  (red  arrows)  separated  by  40  micron  oil  gaps  (green 
arrows)  in  a  non-hydrophobized  capillary.  Fig.  lc  shows  a 
space-time  plot  of  the  BZ  drops  in  Fig.  lb  with  initial  [MA]  = 
40  mM.  With  these  initial  conditions,  the  drops  spontaneously 
oscillate.  All  the  drops  are  initially  synchronized  to  have  the 
same  phase  by  first  suppressing  the  oscillations  via  strong 
illumination  for  several  minutes  and  then  removing  the  light 
from  all  the  drops  simultaneously.  During  the  first  two  oscilla¬ 
tions,  drops  #1  through  #6  oscillate  in-phase  while  drop  #7  is 


out  of  phase.  With  time,  the  remaining  drops  sequentially  fall 
out  of  phase  starting  with  drop  #6,  then  #5,  etc.  By  the  12th 
oscillation,  all  seven  drops  are  out  of  phase  with  each  other. 
However,  at  the  14th  oscillation,  drops  #5  and  #6  synchronize 
in-phase  and  at  the  15th  oscillation,  drops  #4,  #5,  and  #6  are 
synchronized  in-phase.  At  this  time,  drops  #1,  #2,  and  #3  stop 
oscillating  altogether.  During  the  out-of-phase  regime  of  the 
first  15  oscillations,  the  phase-shift  between  adjacent  droplets 
is  never  n  (180°).  For  instance,  the  phase-shift  between  droplets 
#6  and  #7  is  less  than  0.7tt  during  most  of  the  experiment.  As 
previously  reported,  such  a  condition  during  the  out-of-phase 
regime  can  be  considered  as  a  mark  of  weak  inhibitory  coupling.1 
Our  explanation  for  the  evolution  over  time  from  in-phase  to  out-of- 
phase  behavior  for  the  first  15  oscillations  is  that  the  initial  in-phase 
condition,  set  by  the  imposition  of  light,  is  unstable,  and  the  system 
evolves  towards  the  stable  anti-phase  attractor.  However,  over  time 
[MA]  decreases  as  it  is  consumed  in  the  closed  system.  Our 
numerical  studies  of  a  pair  of  interacting  drops  in  an  open  system 
reveal  that  at  low  [MA]  the  diffusive  coupling  of  Br2,  HBr02,  and 
Br02*  leads  to  bistability  between  in-phase  and  out-of-phase  attrac¬ 
tors,  which  explains  why  after  15  oscillations  some  of  the  drops 
adopt  an  in-phase  attractor.2  Finally,  at  lower  [MA]  some  drops  stop 
oscillating.  We  identify  those  drops  with  stationary  Turing  states, 
triggered  by  strong  inhibitory  coupling,  which  are  predicted  to  be 
bistable  with  in-phase  attractors  at  very  low  [MA].2 

Fig.  Id  shows  the  period  as  a  function  of  time  for  each  of  the 
seven  droplets  as  numbered  on  the  space-time  plot.  The  periods 
for  the  droplets  that  oscillate  in-phase  and  for  the  droplets 
oscillating  out-of-phase  with  adjacent  neighbors  evolve  along 
two  different  lines.  At  the  beginning  of  the  experiment,  only  one 
droplet  (#7)  is  out  of  phase  with  the  other  drops  and  lies  on  the 
out-of-phase  line,  which  has  a  longer  period  than  the  in-phase 
drops  at  the  corresponding  time.  Once  a  droplet  moves  from  the 
initially  unstable  in-phase  behavior  to  the  more  stable  out-of- 
phase  attractor,  the  period  increases  and  the  drop  switches  to 
the  upper  line.  This  process  occurs  more  gradually  in  droplet  #1, 
which  does  not  go  back  to  in-phase  interaction.  At  a  later  stage 
some  droplets,  e.g.  #5  and  #6,  return  to  their  earlier  in-phase 
behavior  and  consequently,  their  period  decreases.  Fig.  le  shows 
the  temporal  behavior  of  a  single  drop  that  is  extracted  from  an 
array  of  interacting  drops  (see  Fig.  10a)  with  initial  conditions 
[MA]  =  60  mM,  [Br03_]  =  600  mM  and  [H2S04]  =  160  mM.  The 
intensity  of  transmitted  light,  which  is  proportional  to  the 
ferroin  concentration,  is  also  plotted  as  a  function  of  time.  There 
is  a  noticeable  lengthening  of  the  period  and  change  of  duty 
cycle,  defined  as  the  fraction  of  each  cycle  spent  in  the  oxidized 
state,  with  time.  Oscillations  with  a  long  duty  cycle  of  the 
oxidation  state  are  observed  in  low  [MA]  solutions,  while  the 
oscillations  have  a  short  duty  cycle  at  high  [MA].2,19  Fig.  If  shows 
how  the  ferroin  concentration,  calculated  in  the  point  oscillator 
model  for  a  pair  of  coupled  drops  in  steady  state,  varies  as  a 
function  of  time  for  different  values  of  [MA].  Matching  the 
waveform  from  the  calculation  with  the  observed  waveforms  in 
Fig.  le  allows  us  to  estimate  how  [MA]  varies  as  a  function  of 
time.  For  this  range  of  conditions  we  find  [MA](£  [s])  «  60  mM  — 
0.03  [mM  s_1]  t  [s].  Fig.  lg  shows  the  numerically  calculated 
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period  (based  on  the  FKN  model)  of  a  single  BZ  droplet  as  a 
function  of  malonic  acid  concentration.  The  period  increases  as 
[MA]  decreases;  thus  all  the  observations  in  Fig.  1  are  consistent 
with  a  decreasing  concentration  of  malonic  acid  in  the  drops 
with  time. 

3.2  Modeling  inhibitory  coupling  strength:  finite  element 
analysis 

Previous  experimental  studies1,3,5  of  a  large  number  of  equally 
spaced  droplets  inside  a  cylindrical  capillary  found  anti-phase 
attractors,  in  which  adjacent  drops  oscillated  180°  out  of  phase 
with  each  other.  This  occurred  for  a  wide  range  of  droplet 
diameters  (a),  length  of  oil  gap  between  droplets  (b)  and  inter¬ 
mediate  [MA].  This  anti-phase  behavior  was  explained  to  be  a  result 
of  inhibitory  coupling  of  the  drops,  which  arises  from  the  diffusion 
of  bromine  (Br2)  between  droplets.  For  clarification,  in  the  aqueous 
BZ  reaction,  inhibition  is  carried  out  by  bromide  (Br_),  but  this 
charged  species  does  not  partition  into  the  oil.  However,  Br2,  which 
is  nonpolar,  partitions  readily  into  the  oil  and  diffuses  from  drop  to 
drop.  Once  inside  the  aqueous  phase,  Br2  brominates  MA,  generat¬ 
ing  Br_;  thus  Br2  acts  as  the  carrier  of  inhibition,  and  not  as  the 
inhibitor  itself:  Br2  +  MA  ->  BrMA  +  Br_  +  H+.  Our  results  suggest 
that  this  bromination  can  be  characterized  by  an  effective  rate 
constant  Jcef£  =  10  [s_1]  +  29  [M_1  s_1]  [MA]. 

(a)  Origin  of  coupling:  diffusive  flux.  In  Fig.  2a,  we  present 
the  results  of  solving  a  reaction-diffusion  finite  element  model 
of  the  FKN  equations.  The  model  accounts  for  the  permeation 
of  the  activators  bromine  dioxide  and  bromous  acid,  as  well  as 
bromine,  which  couples  to  the  inhibitory  reaction,  into  the  oil 
separating  the  BZ  drops.  We  modeled  the  drops  as  spheres 
contained  in  a  cylindrical  capillary  whose  inner  diameter 
matches  the  sphere  diameter.  We  also  modeled  the  spheres 
as  lines  in  1-D  and  disks  in  2-D  and  found  only  minor 
quantitative  differences  between  the  1-D,  2-D  and  3-D  solutions. 
To  investigate  how  two  drops  interact,  we  ran  simulations  on 
three  configurations:  two  drops  of  diameter  D  separated  by  0.1D  to 
represent  conditions  of  strong  coupling,  two  drops  of  diameter  D 
separated  by  10D  to  represent  conditions  of  weak  coupling,  and  a 
single  drop  as  a  reference  for  comparison.  In  Fig.  2  the  drops  are 
200  pm  in  diameter  filled  with  BZ  solution  with  [MA]  =  20  mM  in  a 
2-D,  rectangular,  impermeable  container  filled  with  oil. 

The  observation  that  the  period  of  a  pair  of  180°  out-of-phase 
droplets  is  longer  than  that  of  a  pair  of  drops  with  0°  phase  shift,3 
which  was  described  in  Fig.  Id,  can  be  explained  by  examining 
Fig.  2b,  in  which  the  concentration  of  Br2  in  the  center  of  the  left 
drop  and  in  the  oil  gap  separating  the  two  drops  is  plotted  as  a 
function  of  time.  Bromide  (Br“)  inhibits  the  BZ  oscillation.  Injection 
of  a  pulse  of  Br“  into  a  homogeneous  stirred  BZ  solution  lengthens 
the  period  of  oscillation,  because  additional  time  is  required  for 
removing  it,  which  is  required  in  order  to  start  the  autocatalytic  BZ 
step.20,21  However,  it  is  Br2  and  not  Br“  that  is  transmitted  between 
drops.  Br2  is  rapidly  produced  in  a  drop  when  it  undergoes  the 
oxidation  transition.  If  two  drops  are  out-of-phase,  then  the  pulse  of 
Br2  produced  in  one  drop  (the  transmitter)  diffuses  to  its  neighbor 
(the  receiver),  where  it  is  converted  to  Br“,  thereby  delaying  the 
oscillation  of  the  receiving  drop. 
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Fig.  2  Bromine  coupling,  (a)  Mesh  used  in  a  finite  element  calculation  in  2-D 
of  two  200  |!m  drops  of  BZ  fluid  immersed  in  oil.  The  full  FKN  model  was 
solved  on  each  vertex  with  the  following  conditions:  [MA]  =  20  mM,  20  |im  oil 
gap.  The  boundary  condition  at  the  oil/BZ  interface  was  set  by  the  partition 
coefficient,  P^2  =  [ B r2l oil/ [ B r2l bz  =  2.5,  Pq ro2  =  [Br02]Oj[/[Br02]Bz  =  3,  PHBro2  = 
[HBr02]oii/[HBr02]Bz  =  0.  Here  the  partition  coefficient  for  Br02  was  set  to  an 
unphysical  value  3  to  test  the  effect  of  excitatory  coupling  in  simulation.  The 
same  coupling  strength  can  be  obtained  by  lowering  this  value  to  1  and  using 
Fhbi-o2  «  0.01.  A  no-flux  boundary  condition  was  set  at  the  exterior  boundary. 
The  time  corresponds  to  1778  s  in  Fig.  2b.  Arrows  show  flux  of  Br2  and  color 
shows  [Br2],  with  orange  high  and  blue  low.  (b)  [Br2]  as  a  function  of  time  in 
the  oil  between  the  drops  at  the  midpoint  on  the  symmetry  axis  and  in 
the  center  of  the  leftmost  drop.  Two  separations  of  drops  are  calculated; 
2000  |!m  (10D;  10  times  drop  diameter)  and  20  |im  (0.1D;  0.1  of  drop 
diameter).  A  single  drop  (one  drop)  in  the  larger  geometry  is  also  calculated, 
(c)  [Br2]  as  a  function  of  distance  across  the  simulation  box  of  Fig.  2a  for  two 
different  [MA],  Concentrations  are  plotted  at  1778  s,  just  after  the  left  drop 
undergoes  an  oxidation  transition.  The  range  U)  that  the  [Br2]  penetrates  into 
the  neighboring  drop  is  a  measure  of  the  coupling  strength.  7  is  defined  as  the 
distance  it  takes  for  [Br2]  to  drop  halfway  to  its  minimum.  Decreasing  [MA] 
increases  7  as  shown  in  the  inset. 


The  dashed  blue  curve  (“10D  drop”)  in  Fig.  2b  shows  the 
temporal  variation  of  [Br2]  in  the  center  of  drop  1  in  a  closed 
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container  with  drop  2  located  10  diameters  away.  In  a  second 
calculation  involving  a  container  with  a  single  drop,  the  temporal 
variation  of  [Br2]  is  shown  as  the  dashed  green  curve  (“one  drop”). 
The  location  of  drop  1  and  the  size  of  the  container  were  the  same 
as  in  the  “10 D  drop”  calculation,  but  drop  2  was  removed  and 
replaced  with  oil.  The  temporal  variation  of  [Br2]  for  two  drops 
separated  by  10  diameters  (“10 D  drop”)  is  indistinguishable 
from  the  case  of  one  drop  in  a  sealed  container  of  the  same  size 
(“one  drop”),  which  indicates  that  when  drops  are  10  diameters 
apart  they  are  decoupled. 

When  the  drops  are  brought  close  to  each  other  (a  gap  of 
0.1  drop  diameter)  the  Br2  concentrations  are  significantly 
distorted,  as  shown  in  Fig.  2b  “0.1D  drop”.  Compared  to  the 
case  of  drops  far  apart,  during  half  the  cycle  the  level  of  Br2 
inside  drop  1  is  depressed  and  in  the  other  half  of  the  cycle  the 
Br2  level  in  the  same  drop  is  elevated.  The  explanation  is  that  as 
the  drops  are  out  of  phase  with  each  other,  when  drop  1  under¬ 
goes  the  oxidation  transition  at  t-  1778  s,  it  emits  a  large  amount 
of  Br2  at  the  moment  when  the  Br2  level  in  its  neighbor,  drop  2,  is 
low.  Thus  drop  2  acts  as  a  sink  of  Br2  leading  to  a  reduction  in 
[Br2]  in  drop  1  compared  to  the  case  when  the  drops  are  greatly 
separated.  Likewise,  when  drop  2  oxidizes  at  t  =  2400  s,  it  emits  a 
pulse  of  Br2,  raising  the  Br2  level  of  drop  1  above  its  value  when 
the  drops  are  far  apart  and  causing  the  period  of  oscillation  to 
lengthen.  Note  that  the  diffusive  coupling  is  fast  when  the  drops 
are  separated  by  an  oil  gap  of  20  pm;  thus  drop  1  receives  the 
pulse  of  Br2  less  than  1  s  after  it  is  emitted  by  drop  2  at  2400  s. 

Fig.  2a  shows  the  flux  and  concentration  of  Br2  and  Fig.  2b  plots 
the  concentration  as  a  function  of  time.  We  also  measured  the  total 
flux  in  and  out  of  a  drop  as  a  function  of  partition  coefficient  in 
steady  state  conditions  (not  shown).  Under  these  conditions  the 
total  flux  is  zero.  We  found  that  the  flux  out  of  a  drop  increased 
monotonically  with  partition  coefficient  in  steady  state,  however  it 
is  not  clear  how  to  quantitatively  relate  the  chemical  flux  to  the 
degree  of  dynamical  coupling  between  drops. 

(b)  Dimensional  analysis  and  S  parameter.  One  heuristic 
measure  of  the  coupling  strength  is  the  dimensionless  number 
S  =  (PBTD)/(a(a  +  b)Jce ff)  =  /i//ceff,  which  is  the  ratio  of  two  rates: 
p  =  ( PBT  D)/{a(a  +  b )),  the  rate  of  diffusive  transport  between  BZ 
drops  of  the  inhibitor  bromine  separated  by  an  oil  gap,  and  keffl 
the  effective  rate  constant  characterizing  the  consumption  of 
Br2,  which  as  described  previously,  occurs  via  bromination  of 
malonic  acid.  Here  PBr2  is  the  partition  coefficient,  D  is  the 
diffusion  coefficient  of  bromine,  a  is  the  BZ  droplet  size,  b  is  the 
oil  gap  size.  The  derivation  of  p  (assuming  the  drop-oil  inter¬ 
face  is  flat)  is  presented  elsewhere;1  an  alternative  expression 
for  fi  using  the  Derjaguin  approximation  considering  the  actual 
curvature  of  the  drops  is  derived  in  another  work.8 

In  order  to  estimate  Jce ff,  we  plot  the  concentration  of  Br2  as  a 
function  of  distance  in  Fig.  2c  at  t  =  1778  s,  just  after  the  drop  on 
the  left  is  oxidized,  as  depicted  in  Fig.  2a.  If  the  time  for  bromine  to 
diffuse  across  a  drop  is  short  compared  to  the  oscillation  period, 
then  dc(x,t)/dt  can  be  neglected.  For  D  =  10-5  [cm2  s^1]  and 
a  =  1CT2  [cm]  the  time  for  Br2  to  diffuse  across  the  drop,  td  =  10  s, 
is  much  shorter  than  the  oscillation  period,  which  is  greater  than 
300  s.  The  solution  to  the  time-independent  reaction-diffusion 


equation  is  an  exponential  with  a  characteristic  decay  length 
7  =  y/D/keff.  We  define  7  as  the  distance  from  the  water/oil 
interface  to  the  point  that  the  Br2  concentration  is  reduced 
halfway  to  its  minimum.  We  find  7  is  constant  in  time, 
supporting  the  contention  that  dc(x,t)/dt  is  negligible. 
In  Fig.  2c  we  plot  7-2  vs.  malonic  acid  concentration  and 
numerically  find  that  Jcef£  is  linear  in  malonic  acid  concen¬ 
tration;  ke ff  =  (70  M-1  s_1)  ([MA]).  This  result  is  roughly 
consistent  with  the  7-variable  mechanism,  which  suggests  that 
Jce ff  should  be  equal  to  k6  +  k7  =  (10  +  29  [MA]  (M-1))  s_1. 

Two  factors,  one  chemical  and  one  geometrical,  control  how 
strongly  the  bromine  emitted  from  one  drop  influences  the  oscilla¬ 
tion  of  the  receiving  drop.  In  the  limit  of  S  =  ( PD/a(a  +  b))/Jceff  «  1, 
the  drops  are  weakly  coupled.  Chemically,  S  decreases  when  keff  is 
increased.  This  is  accomplished  by  increasing  (MA),  which 
increases  Br2  consumption  inside  a  drop,  leaving  less  Br2  to  diffuse 
between  drops.  Furthermore,  at  high  [MA],  the  Br2  that  is  emitted 
from  one  drop  and  transported  through  the  oil  is  rapidly  consumed 
in  the  receiving  drop  upon  arrival  and  therefore  only  slightly 
increases  the  Br2  concentration  in  the  receiving  drop.  Geometri¬ 
cally,  increasing  the  separation  of  the  drops  weakens  coupling.  As 
transport  is  diffusive,  this  both  increases  the  transport  time  and 
broadens  the  transmitted  pulse,  thereby  reducing  the  diffusive 
transport  between  drops.  In  the  limit  of  S  »  1  the  drops  are 
strongly  coupled.  Decreasing  (MA)  strengthens  coupling  because 
the  reaction  rate  Jce ff  decreases.  When  Jceff  is  small  and  S  is  large, 
bromine  diffuses  from  the  emitting  drop  to  the  receiving  drop 
faster  than  the  receiving  drop  can  consume  the  transmitted 
bromine.  Consequently,  the  concentration  of  bromine  in  the 
receiving  drop  will  be  significantly  increased,  resulting  in  either  a 
phase  shift  in  the  limit  cycle,  or,  for  very  strong  coupling,  distorting 
the  shape  of  the  limit  cycle.  Alternatively,  the  coupling  strength  can 
be  increased  by  decreasing  the  drop  separation,  which  increases 
the  inter-drop  diffusive  transport.  We  previously  observed  that  as 
either  (MA)  or  drop  size  decreases,  more  drops  stop  oscillating  out 
of  phase  and  enter  a  non-oscillatory  Turing  state,  in  which  drops 
are  maintained  far  from  equilibrium  in  either  the  reduced  or 
oxidized  state.1  This  observation  is  consistent  with  theoretical  work 
on  coupled  BZ  oscillators,  which  predicts  that  a  stationary  Turing 
state  is  reached  as  the  inhibitory  coupling  is  increased.5,22 

3.3  Phase  model  for  weakly  coupled  oscillators 

We  now  consider  a  second  measure  of  the  coupling  strength. 
Drops  that  beat  in  synchrony  are  coupled.  Experimentally  and 
numerically  we  have  found  that  for  a  pair  of  drops  there  are  three 
limiting  synchronous  behaviors:  in-phase,  in  which  the  oscillators 
have  the  same  phase,  anti-phase,  in  which  the  oscillators  beat 
180°  out  of  phase13  and  multistable,  in  which  multiple  attractors 
are  possible,  depending  on  initial  conditions.  A  natural  measure 
of  the  strength  of  the  coupling  is  the  degree  to  which  one  drop 
affects  the  phase  of  its  neighbor,  characterized  by  the  phase 
coupling  function,23-25  H.  The  phase  model  is 

Oi  =  ojq  +  H{62  —  $i)  (la) 

02  =  co o  +  H(81  -  e2)  (lb) 
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with  6i  the  phase  of  the  oscillator,  co0  the  frequency  of  each 
oscillator  when  uncoupled,  and  H  the  phase  coupling  function. 
If  H  does  not  change  sign,  then  there  is  only  one  attractor, 
either  in-phase  or  anti-phase,  set  by  the  sign  of  the  derivative  of  H. 
The  number  of  zero  crossings  of  H  sets  the  number  of  multistable 
attractors  (see  Kuramoto,13  eqn  (5.2.18)). 

(a)  Single  species  effect  on  coupling.  Before  we  calculate  the 
phase  coupling  function,  we  first  numerically  explore  the 
qualitative  effect  transmitted  chemicals  have  on  the  synchro¬ 
nization  of  coupled  drops.  We  do  this  by  allowing  only  one 
species  at  a  time  to  diffuse  between  drops  and  numerically 
solve  our  point  oscillator  and  finite  element  models.  The  naive 
expectation  is  that  coupling  through  inhibitory  species  would 
lead  to  out-of-phase  synchronization  and  excitatory  coupling  leads 
to  in-phase  synchrony.  If  a  pair  of  oscillators  are  inhibitory  coupled, 
then  each  delays  the  other,  and  their  phase  difference  will  grow  to 
the  largest  possible  value,  180°.  Such  a  coupling  is  called  “phase 
repulsive”.13  By  the  same  logic,  coupling  of  excitatory  species  would 
lead  to  in-phase  coupling,  because  each  oscillator  stimulates 
the  other,  shortening  their  phase  difference,  referred  to  as 
“phase  attractive”  coupling.13 

We  consider,  one  at  a  time,  transport  of  the  inhibitory 
species,  Br2  and  Br-,  and  the  excitatory  species,  HBr02  and 
Br02#,  in  the  point  oscillator  and  finite  element  models.  In  the 
point  oscillator  model,  if  we  allow  only  Br2  to  transport  between 
drops  and  limit  the  partition  coefficient  PBr2  <  400,  the  drops 
synchronize  out-of-phase.  While  the  measured  value  is  Pby2  ~  2.5, 
for  the  sake  of  completeness  we  consider  Pbx2  >  400,  in  which 
case  a  Turing/in-phase  (depending  on  initial  condition)  mixed 
state  is  found.  For  the  finite  element  model,  the  result  is  qualita¬ 
tively  the  same,  however  the  partition  constant  at  which  the 
behavior  switches  from  out-of-phase  to  in-phase  occurs  at  an  8- 
fold  lower  value,  PBr2  =  50,  than  in  the  point  oscillator  model.  The 
ion  Br_  is  charged  and  therefore  not  expected  to  partition  into  the 
oil.  However,  again  for  the  sake  of  completeness,  we  consider  how 
inter-drop  transport  of  Br~  influences  synchronization.  We  find 
that  in  the  point  oscillator  model,  for  a  partition  coefficient  PBr-2  < 
0.1,  the  drops  synchronize  out-of-phase.  For  larger  partition 
coefficients  the  behavior  first  becomes  bistable  between  out-of- 
phase  and  in-phase  and  then  solely  in-phase.  The  same  qualitative 
behavior  is  observed  for  the  finite  element  model.  These  two 
examples  demonstrate  that  purely  inhibitory  coupling  can  lead  to 
in-phase  synchronization,  contrary  to  our  naive  expectation,  but 
only  when  the  partition  coefficients  are  unphysically  high.  When 
in-phase  synchrony  occurs,  the  drops  are  resetting,  by  which  we 
mean  that  the  oxidation  transition  of  one  drop  induces  an 
immediate  oxidation  transition  in  its  neighbor.  We  note  that  it 
is  possible  in  principle  for  in-phase  synchrony  in  coupled  non¬ 
linear  oscillators  to  occur  without  resetting,  but  we  never  find  this 
in  our  experiments  or  simulations.  The  intuitive  description  is 
that  at  high  inhibitory  coupling  strength,  two  drops  behave  as 
one.  With  only  excitatory  coupling  (restricting  transport  between 
drops  to  either  of  the  two  activators,  HBr02  or  Br02#),  only 
in-phase  synchronization  is  observed. 

In  other  words,  weak  inhibitory  coupling  is  phase  repulsive 
and  excitatory  coupling  is  phase  attractive.  Strong  coupling, 


whether  inhibitory  or  excitatory,  produces  in-phase  synchrony. 
The  dynamical  behavior  of  coupled  drops  is  summarized  in 
Fig.  3a.  For  pure  Br2  coupling,  we  identify  the  boundary  between 
weak  and  strong  coupling  with  the  phase  transition  between  out- 
of-phase  and  in-phase  coupling.  Using  the  values  for  the  point 
oscillator  model  at  the  transition,  a  =  b  =  150  pm,  [MA]  =  200  mM, 
Pbx2  =  400  and  /ceff  =  70  [MA  (M)]  (M_1  s_1),  the  coupling  strength 
S  =  (PD)/(a(a  +  b)Jce ff)  =  0.63;  thus  for  this  case  the  transition  from 
weak  to  strong  coupling  occurs  at  S  «  0.6. 

(b)  Dynamical  phase  boundary  and  S  contour.  To  test  further 
the  conjecture  that  S  is  a  measure  of  coupling  strength,  we 
calculated  the  phase  behavior  using  COMSOL  for  two  drops 
that  are  coupled  solely  by  the  diffusion  of  bromine  through  the  oil. 
The  result  is  plotted  as  a  function  of  [MA]  and  drop  size  for  the  case 
of  oil  gap  equal  to  drop  diameter  and  PBr2  =  2.5,  shown  in  Fig.  3b. 
Superimposed  are  green  contour  lines  of  constant  coupling 
strength,  S.  If  S  is  a  measure  of  coupling  strength,  then  the 
boundary  between  the  anti-phase  (AP)  oscillatory  and  stationary 
Turing  states  should  occur  for  a  constant  value  of  S.  Fig.  3b  shows 
that  the  boundary  corresponds  roughly  to  S  «  0.1,  which  supports 
the  contention  that  the  dimensionless  parameter  S ,  governs  the 
dynamical  behavior  of  the  coupled  BZ  drops. 

Additionally,  for  the  case  of  bromine-only  coupling,  we 
calculated  the  dynamical  phase  boundary  between  anti-phase 
and  stationary  Turing  states  as  a  function  of  drop  size  a  and  oil 
gap  size  b  using  COMSOL  and  superimposed  S  contours  as 
shown  in  Fig.  3c.  While  the  construction  of  the  point  model 
requires  that  the  phase  boundary  is  a  function  of  a  x  (a  +  b),  the 
functional  form  is  not  specified.  In  particular,  there  is  no 
reason  to  suppose  that  the  phase  behavior  scales  as  inversely 
proportional  to  a  x  (a  +  b),  as  is  demonstrated  in  Fig.  3c. 
Furthermore,  there  is  no  such  constraint  on  the  phase  boundary  in 
the  finite  element  models,  as  those  equations  do  not  have  a  unique 
form  of  non-dimensionalization.  We  interpret  the  results  of  Fig.  3c 
to  mean  that  small  drops  are  strongly  coupled  and  that  drop-drop 
coupling  is  relatively  insensitive  to  drop  separation.  Experimental 
study  of  the  dynamical  behavior  of  drops  as  a  function  of  size  and 
separation  will  be  the  subject  of  future  work. 

(c)  Phase  coupling  function  H.  Next  we  calculate  H,  the 
phase  coupling  function.  Using  the  point  oscillator  model,  we 
numerically  calculate  the  phase  coupling  function  (see  eqn  (5.2.17b) 
of  Kuramoto,13  or  eqn  (10.15)  of  Izhikevich25)  for  two  identical 
coupled  BZ  drops  as  a  function  of  drop  diameter,  drop  separation, 
and  [MA].  We  consider  that  the  following  3  components  of  the 
Turing  point  oscillator  model  have  non-zero  partition  coefficients; 
PBr2  =  2.5,  Pbvo2  =  1,  and  PHbkd2  =  0.1.  The  first  two  were  estimated 
from  experiment;  the  last  was  obtained  by  varying  PHbkd2  until  the 
best  fit  between  experimental  and  calculated  dynamical  phase 
diagrams  was  obtained.  Thus  PHbio2  =  0.1  is  a  prediction  of  our 
point  oscillator  model.  Using  the  finite  element  model  we  obtained 
^Wo2=  0.01.  While  there  were  quantitative  differences  between 
the  two  models,  their  qualitative  behavior  was  similar. 

There  are  two  ways  to  calculate  H.  One  is  indirectly,  as  done 
by  Kuramoto,13,24  using  the  phase  response  curve  of  a  single 
drop,  which  in  principle  could  be  calculated  from  the  point 
oscillator  model.  Instead  we  obtain  H  directly  by  numerically 
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Fig.  3  Dynamical  phase  boundary  and  coupling  strength  S.  (a)  Calculated  in  the  point-oscillator  model  with  a  =  b  =  150  pirn,  and  [MA]  =  200  mM.  Only 
one  chemical  species  at  a  time  is  allowed  to  diffuse  between  two  drops.  Four  species  are  considered;  two  inhibitory  and  two  excitatory.  Phase  behavior; 
AP;  anti-phase  oscillation.  Turing  (stationary);  non-steady  state.  IP;  in-phase  oscillation.  "AP/IP"  indicates  a  multistable  state.  The  results  are  qualitatively 
the  same  as  those  obtained  from  the  finite  element  analysis  (not  shown).  The  partition  coefficient  at  the  phase  transition  is  indicated  for  the  point 
oscillator  model.  For  Br2  in  the  point  oscillator  model,  the  coupling  strength  S  =  0.6  at  the  phase  transition  point,  (b)  Phase  diagram  as  a  function  of 
drop  size  and  malonic  acid  concentration  (log  scale)  calculated  using  COMSOL  (1-D)  for  two  BZ  drops  coupled  through  only  bromine  with  partition 
coefficient  2.5  and  a  =  b.  The  dark  red  region  is  in  a  stationary  Turing  state,  while  the  blue  region  is  anti-phase  oscillatory.  Superimposed  green  lines 
are  contour  lines  of  coupling  strength  S;  the  numbers  are  the  values  of  S.  (c)  Similar  to  Fig.  3b,  phase  diagram  as  a  function  of  drop  size  and  oil  separation 
with  [MA]  =  60  mM. 


calculating  the  phase  shift  that  one  drop  imposes  on  its 
neighbor26  after  one  period.  Consider  two  point  oscillators  that 
are  uncoupled  (H  =  0)  and  have  a  constant  phase  shift  6.  Since 
they  are  uncoupled  6t  =  co0.  The  ferroin  concentration  of  these 
drops  is  plotted  as  a  function  of  time  in  Fig.  4a  (blue  curves).  At 
an  arbitrary  time,  labeled  t  =  0  in  Fig.  4a,  we  couple  the  two 
drops  via  the  oil  and  plot  the  ferroin  concentration  for  the 
coupled  drops  as  a  function  of  time  (red  curves).  After  one 
period  of  the  uncoupled  drop  (t),  the  phase  shift  of  drop  1  (AO), 
relative  to  its  phase  when  uncoupled  is  calculated  and  in  Fig.  4b 
is  plotted  as  a  function  of  0  for  four  different  values  of  malonic 
acid.  A6/t  is  equal  to  the  phase  coupling  function  called  r  by 
Kuramoto  in  eqn  (5.2.17)24  and  //by  Izhikevich  in  eqn  (10. 16). 25 
We  use  units  of  phase  that  vary  between  0  and  1.  The  phase 
shift  is  unambiguous,  because  the  coupled  and  uncoupled 
drops  move  on  the  same  limit  cycle,  x(t),  where  x  represents 
the  7  chemical  variables  in  the  FKN  model.  We  numerically 
verified  that  A 0  is  independent  of  the  time  when  the  two  drops 


are  coupled;  for  example,  if  the  drops  in  Fig.  4a  are  coupled 
beginning  at  t  =  50  s  or  t  =  —50  s,  A 6  is  the  same  as  calculated 
for  t  =  0.  In  Fig.  4a,  the  time  when  A 6  is  calculated  is  indicated 
as  a  dashed  line.  The  value  of  A 6  is  slightly  less  than  AO',  the 
phase  difference  between  the  expected  oxidation  transition  for 
an  uncoupled  drop  and  for  the  same  drop  after  it  is  coupled. 
However,  AO'  is  notable  because  it  is  readily  measurable.26  As 
shown  in  Fig.  1,  we  measure  the  ferroin  concentration  during 
the  entire  cycle,  from  which  it  is  possible  in  principle  to  extract 
A 6'.  However,  in  practice,  while  numerically  calculating  A 6  is 
easy,  experimentally  it  is  subject  to  a  large  uncertainty.  The  sign 
of  the  phase  shift  tells  whether  or  not  the  coupling  is  phase 
attractive  (negative)  or  phase  repulsive  (positive).  As  a  measure 
of  the  coupling  strength,  we  consider  the  value  of  A 6  when  0  = 
0.5,  that  is  the  phase  shift,  AO,  induced  by  coupling  two  drops 
whose  initial  phase  difference  is  0  =  0.5,  or  180°.  This  particular 
choice  for  6  is  somewhat  arbitraiy.  Using  this  definition,  we 
find  that  the  repulsive  coupling  strength  increases  as  the  drop 
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Fig.  4  Calculated  coupling  strength,  (a)  Ferroin  concentration  \/s.  time 
calculated  using  the  FKN  point-oscillator  model  for  two  drops  such  as 
illustrated  in  Fig.  2a  in  the  out-of-phase  condition.  Thin  blue  lines: 
decoupled  drops.  Thick  red  lines:  coupled  drops.  Drops  are  coupled  at 
t  =  0  and  the  phase  shift  [AO)  between  the  coupled  and  decoupled  drops  is 
calculated  after  one  period  (t)  of  the  uncoupled  system  as  a  function  of  the 
initial  phase  shift  {0)  between  the  decoupled  drops.  Phase  is  normalized  to 
1.  (b)  Phase  shift  after  one  period  of  coupling  [AO)  vs.  initial  phase  shift  [0) 
for  different  values  of  [MA]  using  the  point  oscillator  model,  with  PHBro2  = 
0.1,  PBr2  =  2.5,  PB ro2  =  1,  a  =  b  =  150  m.  A 6  increases  with  0  and  increases 
as  [MA]  decreases,  (c)  Phase  diagram  of  attractor  behavior  as  function  of 
drop  size  and  [MA]  with  the  same  partition  coefficients  as  in  Fig.  4b.  The 
inset  contour  plot  shows  how  the  magnitude  of  AO  for  0  =  n  varies  as  a 
function  of  [MA]  and  drop  size.  IP:  in-phase  synchronization.  AP/IP: 
bistability.  T:  stationary  Turing  state,  in  which  some  drops  stop  oscillating, 
(d)  Dynamics  of  the  phase  difference  between  two  oscillators  (A 0)  calcu¬ 
lated  with  the  point  oscillator  (dotted  curves)  and  phase  (solid  curves) 
models  using  [MA]  =  600  mM  and  the  same  geometry  and  partition 
coefficients  as  in  Fig.  4b.  A (f)t=0  =  0.1,  0.3,  0.5  for  the  blue,  red,  and  green 
curves,  respectively.  Phase  varies  between  0  and  1. 


size  decreases,  as  the  drop  separation  decreases,  and  as  the 
malonic  acid  concentration  decreases,  consistent  with  our  previous 
measure  of  coupling  strength,  S.  We  find  that  the  coupling  is  purely 
repulsive  for  a  wide  range  of  conditions:  200  mM  <  [MA]  <  2  M 
and  drop  size  >100  pm,  consistent  with  our  previous  observations 
with  such  drops.1  For  smaller  drops  or  low  [MA],  bistable  attractors, 
in-phase  attractors  and  stationary  Turing  states  appear.  In  Fig.  4b, 
for  the  case  [MA]  =  20  mM,  the  coupling  function  crosses  zero, 
rendering  it  bistable,  with  an  in-phase  attractor  above  6  >  0.9 
and  anti-phase  otherwise.  The  dynamical  phase  behavior  is 
summarized  in  Fig.  4c. 

(d)  Validity  of  phase  model.  We  also  studied  the  dynamics  of 
phase  coupling  for  a  pair  of  droplets  using  the  point  oscillator 
model  (Fig.  4d).  Similar  to  eqn  (10.13)  to  (10.17)  in  Izhikevich27 
and  Section  5.2.2  in  Kuramoto,13  we  define  0$)  =  tfr  + 
i-  1,  2  for  phases  of  the  two  droplets,  with  the  first  term  the 
phase  of  free-running  oscillation  (t  is  the  free-running  period) 
and  the  second  term  the  slow  phase  deviation  induced  by 
coupling  these  two  droplets.  In  Fig.  4d  we  plot  A cj)(t)  =  ^{t )  — 
(j>2{t)  as  calculated  from  the  point  oscillator  model  and  the 
phase  model  (eqn  (la)  and  (b))  for  various  initial  conditions, 
A(j)t=0  -  0.1,  0.3,  0.5.  For  both  models  the  phase  difference 
evolves  from  the  initial  condition  to  a  steady  phase  deviation, 
oo  =0.5,  corresponding  to  anti-phase  coupling.  The  results 
for  the  point  oscillator  model  give  the  instantaneous  phase 
difference  (dotted  lines),  while  the  phase  oscillator  model  (solid 
lines)  represents  only  the  slow  dynamics,  as  the  coupling  function 
H  was  obtained  by  averaging  the  coupling  between  drops  over  one 
period.  The  fact  that  A </>  corresponds  for  the  point  and  phase 
models  establishes  the  validity  of  the  phase  model  for  BZ  droplets 
in  the  weak  coupling  limit. 

3.4  Excitatoiy  coupling  and  in-phase  behaviors 

(a)  Trigger  wave.  Using  [MA]  =  40  mM  and  droplets  in  contact, 
or  using  [MA]  =  20  mM  and  more  closely  spaced  droplets  than  in 
Fig.  1,  we  were  able  to  obtain  a  stable  in-phase  behavior  of  large 
arrays  of  ID  droplets  as  observed  in  Fig.  5(a-d).  We  were  unable  to 
observe  such  patterns  for  40  mM  <  [MA]  <  2  M.  In  contrast 
with  previously  reported  unstable  in-phase  patterns,3  the  stable 
patterns  found  here  often  (but  not  always)  initially  begin  with  zero 
phase  shift  between  adjacent  drops  and  subsequently  develop 
a  small,  but  constant  phase  shift  between  adjacent  droplets, 
creating  a  propagating  BZ  wave,  analogous  to  the  behavior 
observed  in  a  continuous  aqueous  BZ  solution. 

We  simulated  the  observed  wave  pattern  using  COMSOL  to 
model  a  chain  of  2D  disks  in  oil.  We  postulated  that  the  waves 
observed  experimentally  in  Fig.  5b  and  d  were  trigger  waves.  To 
numerically  explore  this  possibility  we  set  one  drop  on  the  edge 
of  the  simulation  box  to  have  a  higher  concentration  of  sulfuric 
acid  than  the  other  drops,  causing  that  drop  to  oscillate  faster 
than  the  other  drops  and  thus  act  as  a  trigger  (Fig.  5g  and  h). 
We  compared  experiments  showing  traveling  waves  in  drops 
in  cylindrical  capillaries  with  simulation  and  found  similar 
behavior.  To  validate  the  simulations,  we  modeled  a  continuous 
BZ  solution  confined  in  the  same  capillary  and  compared  with 
experiment.  Again  the  simulation  and  experiment  were  similar, 
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Fig.  5  In-phase  patterns:  trigger  wave  experiments  \/s.  simulations,  (a,  c) 
photographs  of  drops  in  capillary,  (b,  d)  Space-time  plots,  (a,  b)  [MA]  = 
20  mM,  (c,  d)  [MA]  =  40  mM.  Droplet  diameter  is  (a,  b)  150  pim  and  (c,  d) 
97  nm;  distance  between  droplets  is  (a,  b)  36  jam  and  (c,  d)  touching 
droplets,  (e)  Experiment:  single  phase  BZ  solution  ([MA]  =  20  mM).  Trigger 
wave  speed  2  mm  min-1,  (f)  Experiment:  BZ  droplets  (100  [im),  trigger 
wave  speed  1  mm  min-1,  (g)  Simulated  in  COMSOL,  single  phase  using 
[MA]  =  20  mM.  Trigger  wave  speed,  4  mm  min-1,  (h)  Simulated  in  COMSOL, 
BZ  droplets  (200  |im),  trigger  wave  speed  2  mm  min-1  The  topmost  droplet 
has  slightly  higher  [H+]  =  0.2  M  than  the  rest  of  the  droplets  (0.16  M)  to  act  as 
an  intrinsically  faster  "trigger".  PH Bro2  =  0-  Pbv2  =  2.5,  PBro2  =  5.  Identical 
results  are  obtained  with  PB r2  =  2.5,  PHBro2  «  0.01  and  PBro2  =  1. 


as  illustrated  in  Fig.  5(e-h).  The  one  discrepancy  was  that  the 
frequency  of  the  oscillation  in  the  simulation  was  about  twice 
that  found  in  the  experiment.  This  is  a  shortcoming  of  the  FKN 
model,  and  not  a  consequence  of  the  finite  element  aspect  of 
the  calculation.  This  discrepancy  in  oscillation  frequency  led  to 
the  wave  speed  in  simulation  being  twice  as  fast  as  in  experi¬ 
ment;  however,  the  calculated  wavelength  agreed  with  experi¬ 
ment.  In  both  simulation  and  experiment,  the  speed  of  the 


trigger  wave  in  the  continuous  BZ  fluid  was  twice  the  speed  of  the 
wave  in  the  chain  of  drops.  The  oil  gap  diminishes  the  diffusive  flux 
in  comparison  to  the  continuous  BZ  solution,  offering  a  possible 
explanation  for  the  slower  propagation  speed  in  the  drops. 

To  explore  the  origin  of  the  temporal  evolution  of  the  phase 
between  adjacent  drops  from  initially  zero  phase  to  a  constant 
phase  as  shown  in  Fig.  5,  we  repeated  the  trigger-wave  simula¬ 
tions  of  Fig.  5g  and  h  as  a  function  of  inter-drop  coupling 
strength.  In  Fig.  6a  we  simulate  a  line  of  identical  drops,  except  for 
one  drop  on  the  end,  which  has  a  higher  frequency.  The  drop 
spacing  is  60  microns,  large  enough  that  the  fast  drop  does  not 
launch  a  trigger  wave,  but  the  drops  are  still  sufficiently  coupled  to 
oscillate  in  phase.  Simulated  light  was  shone  on  the  middle  drop 
with  sufficient  intensity  to  suppress  oscillation  in  that  drop.  One 
observes  that  drops  immediately  adjacent  are  slightly  influenced 
by  the  increased  bromine  produced  by  the  illuminated  drop. 


experiment 


‘''time 


Fig.  6  Traveling  waves  and  coupling  strength.  Blue  arrows  indicate  light- 
suppressed  BZ  drops,  (a,  b)  COMSOL  simulations  ([MA]  =  20  mM,  PHBro2  =  0- 
PB r2  =  2.5,  PBr q2  =  5):  (a)  drop  size  a  =  160  jam,  oil  gap  b  =  60  (jm.  (b)  a  = 
200  urn,  b  =  20  pirn,  light  intensity  0.01  s-1  (pseudo  concentration  in  the 
model).  In  both  (a)  and  (b)  the  first  drop  (top,  [H+]  =  0.22  M)  is  made  to 
oscillate  faster  than  the  other  drops  ([H+]  =  0.16  M).  The  light-suppressed 
drop  in  the  middle  separates  the  drop  array  into  two  parts  in  each  case. 
Trigger  wave  behavior  is  observed  only  in  (b),  where  the  inter-drop  coupling 
is  stronger,  but  in-phase  attractors  are  observed  for  both  (a)  and  (b). 
(c)  Experiment  ([MA]  =  20  mM,  a  =  100  |im,  b  =  10  |im).  Although  the  alb 
ratio  in  this  experiment  is  the  same  as  in  Fig.  6b,  the  result  resembles  Fig.  6a. 
The  experiment  started  with  an  in-phase  attractor  with  zero  phase  shift 
between  drops  and  evolved  into  a  trigger  wave. 
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Fig.  7  Space-time  plots  of  mixed  mode  patterns,  (a)  Droplets  are  170  mm 
in  diameter  separated  by  70  \im,  [MA]  =  20  mM.  (b)  Droplets  of  100  |im 
diameter  are  in  contact,  [MA]  =  40  mM.  (c)  Droplets  of  70  pirn  diameter  are 
in  contact,  [MA]  =  60  mM.  Before  starting  each  experiment,  the  droplets 
were  synchronized  in-phase  by  illuminating  the  capillary  with  intense 
450  nm  light  for  about  10  min.  In  the  space-time  plots,  long  double 
arrows  mark  periods  of  at  least  two  in-phase  neighbor  droplets;  short 
arrows,  half  the  size  of  the  long  arrows,  mark  the  period  of  a  single  droplet. 


In  Fig.  6b  the  inter-drop  separation  is  reduced  to  20  pm.  Now  a 
trigger  wave  is  formed,  but  does  not  propagate  across  the  light- 
induced  gap.  From  these  simulations  we  conclude  that  trigger 
waves  require  a  stronger  coupling  between  drops  than  do  in-phase 
attractors  with  zero  phase  shift  between  neighbors.  Fig.  6c  shows 
an  equivalent  experiment.  The  drop  in  the  center  of  the  capillary 
is  illuminated,  suppressing  its  oscillation.  The  drops  initially 
oscillate  with  zero  phase  shift  between  neighbors,  but  with  time, 
a  traveling  wave  emerges  as  the  malonic  acid  concentration 
decreases  and  the  coupling  strength  consequently  increases. 

(b)  Cluster  patterns.  In  Fig.  7a,  we  present  a  space-time  plot 
highlighting  three  pairs  of  in-phase  droplets  numbered  1-6. 
The  pattern  consists  of  a  row  of  drops  in  which  at  a  particular 
instant  in  time  all  drops  form  pairs  where  drops  in  a  pair  have 
the  same  phase,  but  each  pair  is  180°  out-of-phase  with  its 
neighboring  pairs,  i.e.  drops  #1  and  #2  have  0°  phase,  drops  #3 
and  #4  have  180°  phase  and  drops  #5  and  #6  have  0°  phase.  We 
refer  to  this  pattern  as  a  “local  in-phase,  global  out-of-phase” 
attractor.  In  Fig.  7b  we  observe  defects  in  this  synchronization 
pattern  that  occur  when  there  is  an  additional  drop  separating 
two  pairs,  i.e.  drops  #1  and  #2  have  0°  phase,  drops  #4  and  #5  have 
180°  phase  and  drop  #3  is  the  defect.  In  this  case,  the  droplet 
separating  the  two  pairs  (drop  #3)  alternately  synchronizes 
in-phase  with  the  neighboring  pairs  to  each  side;  a  feat  that  is 
only  possible  when  the  middle  droplet  oscillates  with  twice  the 
frequency  of  its  neighbors.  As  in  the  case  of  the  stable  in-phase 
patterns  in  Fig.  5,  the  patterns  observed  in  Fig.  7  require  low  [MA] 
and/or  droplets  with  small  separation.  We  observed  these  local 
in-phase  -  global  out-of-phase  modes  for  [MA]  up  to  60  mM,  as 
observed  in  Fig.  7c.  At  higher  [MA],  even  touching  droplets 
yield  only  out-of-phase  patterns.1  Typically,  as  the  BZ  reaction 
proceeds  towards  equilibrium,  oscillations  eventually  give  way  to  a 
stationary  Turing  pattern  with  some  drops  locked  in  the  oxidized 
state  and  others  in  the  reduced  state. 


Often,  the  local  in-phase  -  global  out-of-phase  pattern  is  mixed 
with  “pure”  out-of-phase  behavior  or  with  stationary  droplets.  In 
Fig.  7a,  on  either  side  of  the  pairs  of  highlighted  arrowed  regions, 
we  observe  drops  where  nearest  neighbors  are  out-of-phase,  while 
in  Fig.  7b  the  right  side  of  the  space-time  plot  shows  stationary 
droplets  while  the  left  side  remains  oscillatory.  When  the  pattern 
includes  stationary  droplets,  as  observed  in  Fig.  7b  and  c,  the 
stationary  droplets  produce  a  variety  of  complex,  often  symmetric, 
patterns.  In  Fig.  7c  we  observe  periodic  stationary  Turing  patterns 
with  wavelengths  of  4  drops;  in  Fig.  7b  the  Turing  wavelength  is 
3  drops.  The  most  commonly  observed  in-phase  clusters  consist 
of  a  pair  of  droplets,  and  the  number  of  droplet  clusters  drops  off 
sharply  with  the  size  of  the  cluster. 

All  our  models  (Vanag  and  Epstein,2  Fig.  4)  predict  that  if 
drops  are  in  a  stationary  Turing  state  and  the  coupling  is 
increased,  then  the  drops  should  resume  oscillating,  but  in-phase. 
We  have  not  observed  this  behavior  for  drops  in  ID  capillaries. 
At  this  time,  it  is  not  clear  if  this  is  because  the  theoiy  is  wrong, 
that  we  have  not  found  the  right  experimental  conditions,  or, 
most  probably,  the  assumption  that  the  closed  system  as  a 
function  of  time  corresponds  to  an  open  system  with  decreasing 
malonic  acid  concentration  fails. 

(c)  Transport  of  BZ  species  through  the  oil.  The  occurrence  of 
stable  in-phase  patterns  suggests  excitatory  coupling  between 
droplets.  The  permeability  and  chemical  reactivity  of  the  BZ 
products  in  the  perfluorinated  oil  that  we  use  have  not  been  fully 
characterized.  This  oil  is  slightly  polar,  with  a  dielectric  constant 
of  5.8,  whereas  octane,  which  was  employed  in  earlier  studies  of 
BZ  emulsions,5  has  a  dielectric  constant  of  1.9.28  As  we  show 
below,  this  perfluorinated  oil  acts  as  a  transport  medium  for 
slightly  polar  molecules  such  as  o^halogens  and  weak  acids.  To 
investigate  the  transport  of  weak  acids,  we  placed  two  water 
droplets  separated  by  a  fluorinated  oil  gap  inside  a  hydrophobized 
capillary.  One  droplet  contained  a  ~  1  M  solution  of  a  weak  acid, 
while  the  other  had  a  6  mM  solution  of  the  pH  indicator  thymol 
blue  at  pH  8.29  We  studied  chloroacetic,  iodoacetic,  and  malonic 
acids.  Equal  volumes  of  the  thymol  blue  solution  and  any  of  these 
acids  gave  an  intense  red-pink  solution.  Fig.  8  shows  a  typical 
example  of  these  experiments  with  malonic  acid.  The  drop 
containing  thymol  blue  (on  the  left)  becomes  darker  as  time  goes 
on,  indicating  the  transport  of  MA  through  the  fluorinated  oil. 
Moreover,  the  oil  also  becomes  darker  in  time.  This  demonstrates 
that  a  small  amount  of  the  indicator  partitions  into  the  oil  phase, 
which  can  be  explained  by  the  fact  that  thymol  blue  is  an 
uncharged  weak  acid  in  its  completely  protonated  red  form. 
Similar  results  were  obtained  for  chloroacetic  and  iodoacetic 
acids.  From  these  observations,  we  conclude  that  mass  exchange 
of  uncharged  weak  acids,  like  MA  and  HBr02,  takes  place  between 
droplets  in  this  fluorinated  oil.  In  contrast  to  the  weak  acids, 
similar  experiments  with  sulfuric  acid,  a  strong  acid,  gave  no 
evidence  of  permeation  through  the  oil. 

To  estimate  the  interdrop  transport  of  bromine  dioxide 
radical  (Br02#),  which  is  too  reactive  to  employ  in  these  experi¬ 
ments,  we  used  the  chemically  similar,  but  stable  molecule, 
chlorine  dioxide,  C102*.  As  described  in  the  Methods  section, 
Cl02#  was  generated  in  drops  by  mixing  sodium  chlorite  with 
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Fig.  8  Activator  permeation.  From  top  to  bottom,  three  pictures  of  the 
same  capillary  (0.5  mm  internal  diameter)  with  increasing  time.  At  the  left,  a 
droplet  of  thymol  blue,  which  darkens  with  time;  at  the  right,  a  reservoir  of 
1.3  M  malonic  acid.  Fluorinated  oil  HFE-7500  in  the  center  also  becomes 
darker  in  time,  revealing  the  presence  of  the  nonionic  form  of  thymol  blue 
in  the  oil.  The  capillary  is  observed  with  510  nm  light. 


nitric  acid.  It  took  about  20  minutes  for  the  [C102#]  concen¬ 
tration  to  reach  a  maximum  inside  a  drop  and  about  50  minutes 
for  the  [C102*]  in  the  oil  to  exceed  that  in  the  drop.  In  contrast, 
the  transport  of  the  inhibitor,  Br2,  from  drop  to  oil  is  immea¬ 
surably  rapid.  In  the  30  seconds  that  it  takes  to  load  the  capillary 
with  droplets,  remove  it  from  the  microfluidic  chip  and  place  it 
in  the  microscope  for  observation,  a  significant  fraction  of  the 
bromine  has  left  the  drop  and  come  to  equilibrium  with  the  oil 
drop.  These  results  suggest  that  the  rate  of  Br02*  transport 
between  droplets  is  two  or  more  orders  of  magnitude  slower 
than  Br2  transport  through  the  same  interface. 

The  results  described  in  the  two  previous  paragraphs  demon¬ 
strate  that  the  transport  of  Br2  from  the  water  through  the 
fluorinated  oil  is  much  faster  than  the  transport  of  all  the  other 
species  tested  here,  i.e.,  Br02#,  MA  and  the  other  weak  acids.  The 
permeation  coefficient,  defined  as  the  product  PD,30  governs  the 
flux  of  chemicals  between  drops,  where  P  is  the  partition  coeffi¬ 
cient  and  D  is  the  diffusion  coefficient.  D  varies  only  slightly 
between  molecular  species  and  cannot  account  for  the  large 
difference  between  transport  rates  of  the  different  species.  How¬ 
ever,  P  can  vary  by  orders  of  magnitude  among  the  BZ  chemical 
species.  For  Br2  we  have  previously  measured  P  in  this  oil  as  2.5, 1 
while  for  C102#  we  measured  the  partition  coefficient  between  oil 
and  water  spectrophotometrically  after  extraction  in  the  bulk, 
yielding  P  =  1.2.  The  value  of  the  permeability  for  Br2  is  at  least 
twice  the  value  of  any  of  the  other  tested  chemical  species. 


The  experiments  described  in  this  section  strongly  support 
our  hypothesis  that  excitatory  coupling  caused  by  the  transport 
of  Br02#  and  HBr02  through  the  oil  plays  a  major  role  in 
facilitating  the  observed  in-phase  behavior  at  low  [MA].  In  the 
next  section  we  describe  experiments  based  on  this  hypothesis 
that  are  designed  to  enhance  excitatory  coupling. 

(d)  Excitatoiy  coupling  enhancement:  observation  and  con¬ 
sequences.  Two  ways  were  explored  to  enhance  excitatory 
coupling.  First,  decreasing  the  concentration  of  Br2  should 
weaken  the  inhibitory  coupling,  thereby  making  the  excitatory 
coupling  more  effective.  Second,  increasing  the  concentration 
of  HBr02  and/or  Br02*  should  directly  strengthen  the  excitatory 
coupling  between  droplets. 

Exploring  the  first  approach,  we  confirmed  that  MA  and  Br2 
are  able  to  react  in  the  fluorinated  oil,  probably  to  give  dibromo- 
malonic  acid  via  a  classic  addition  of  bromine  to  a  double  bond  in 
a  nonaqueous  solvent.31  This  reaction  removes  Br2  directly  and 
thereby  diminishes  the  production  of  bromide,  the  major  BZ 
inhibitory  species,  in  the  water  drop.  In  an  earlier  study,1  we 
reported  out-of-phase  patterns  at  [MA]  =  2.8  M  and  moderately 
spaced  droplets.  By  increasing  [MA]  to  2.9  M  and  placing  the 
droplets  in  contact,  we  were  able  to  generate  in-phase  trigger  wave 
oscillatory  behavior  that  was  maintained  for  about  a  dozen  cycles, 
as  shown  in  Fig.  9.  Neither  our  point  oscillator  nor  our  finite 
element  model  predicts  in-phase  oscillations  at  [MA]  >  2.8  M. 
Our  numerical  models  do  not  account  for  any  chemical  reactions 
in  the  oil  phase.  However,  we  expect  that  eliminating  the  inhibi¬ 
tory  species,  but  not  the  excitatory  species  from  the  oil  phase 
will  lead  to  in-phase  oscillation.  As  a  function  of  time,  the 
[MA]  concentration  decreases,  and  below  a  threshold  of  2.8  M 
the  in-phase  oscillations  will  transform  to  out-of-phase  behavior. 
As  in  Fig.  1,  in-phase  droplets  oscillate  with  a  shorter  period  than 
when  they  oscillate  out-of-phase.  Another  interesting  observation 
is  that  the  speeds  of  the  trigger  waves  associated  with  the  in-phase 
patterns  at  [MA]  =  40  mM  (Fig.  5)  and  2.9  M  (Fig.  9),  are  the  same. 

Next,  we  examine  the  effect  of  direct  enhancement  of  the 
excitatory  coupling  by  increasing  the  activator  concentration  in 
the  oil,  via  the  aqueous  reaction  HBr02  +  Br03_  +  H+  <->  2Br02  + 
H20,  which  is  a  key  step  in  the  FKN  mechanism.10  We  expect 
that  Br02#  is  more  soluble  in  oil  than  the  more  polar  HBr02 
(although  it,  too,  may  permeate  into  the  oil);  therefore  increas¬ 
ing  the  bromate  and  acid  concentrations  should  shift  the 


Fig.  9  High  malonic  acid,  (a)  Droplets  of  65  nm  diameter  are  in  contact, 
[MA]  =  2.9  M.  (b)  Space-time  plot  showing  a  transient  in-phase  trigger 
wave  pattern  followed  by  an  out-of-phase  pattern.  The  moments  when 
the  oxidized  state  flash  occurs  are  enhanced  for  clarity. 
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Fig.  10  Activator  coupling,  (a)  Space-time  plot  for  [MA]  =  60  mM,  [Br03~]  = 
600  mM  and  [H2S04]  =  160  mM.  Droplets  are  170  |im  in  diameter  and 
190  |!m  apart.  The  total  time  of  the  space-time  plot  is  2000  s. 
(b)  Distribution  of  phase-shifts  of  pairs  of  neighboring  droplets  under 
the  conditions  of  (a)  (black  bars),  compared  with  the  distribution  when 
[MA]  =  60  mM,  [Br03“]  =  300  mM  and  [H2S04]  =  80  mM  (shaded  bars). 
Inset  suggests  a  possible  mechanism  of  excitatory  coupling. 

equilibrium  to  the  right,  increasing  [Br02#]  and  thereby 
strengthening  the  excitatory  coupling  through  the  oil. 

Our  earlier  work1  found  that  with  concentrations  of  [H2S04]  = 
80  mM  and  [Br03_]  =  300  mM,  a  60  mM  malonic  acid  concen¬ 
tration  produced  robust  anti-phase  patterns.  Fig.  10a  shows  a 
space-time  plot  at  [MA]  =  60  mM  and  twice  the  previous  concentra¬ 
tions  of  bromate  and  sulfuric  acid.  We  ran  this  experiment  without 
initially  synchronizing  the  droplets  in-phase.  Fig.  10  compares 
the  distribution  of  phase  shifts  between  neighboring  droplets 
under  stationary  conditions1  with  the  two  sets  of  acid  and 
bromate  concentrations.  We  argue  that  the  relative  strength 
of  excitatory  to  inhibitory  coupling  correlates  with  the  phase 
difference  between  coupled  oscillating  drops;  inhibitory-dominated 
coupling  leads  to  phase  repulsion,  in  which  neighbors  have 
phase  differences  of  half  a  cycle,  whereas  excitatory-dominated 


coupling  results  in  phase  attraction,  where  nearest  neighbors 
have  the  same  phase.  The  existence  of  in-phase  pairs  is  clear  at 
the  higher  bromate  and  acidity  conditions;  therefore  increasing 
bromate  increases  excitatory  coupling. 

4.  Conclusions 

While  we  have  not  been  able  to  manipulate  the  inhibitory  and 
excitatory  coupling  completely  independently,  the  present  experi¬ 
mental  system  affords  considerably  more  control  over  coupling 
than  experiments  in  which  interaction  between  oscillators  occurs 
via  an  aqueous  phase,  so  that  all  species  participate  in  proportion 
to  their  concentrations,4  or  systems  in  which  the  inhibitor 
dominates  the  coupling  between  droplets.1,3,5  Numerical  investi¬ 
gations  of  the  complex  patterns  that  we  observe  suggest  that  both 
excitatory  and  inhibitory  coupling  are  involved. 

The  fact  that  our  experimental  system  is  closed  prevents  us 
from  obtaining  true  attractors,  but  the  persistence  of  patterns 
that  exhibit  the  same  qualitative  behavior  for  many  cycles  of 
oscillation  suggests  that  these  behaviors  would  be  stable  if  we 
were  able  to  maintain  constant  conditions  indefinitely.  This 
conjecture  is  supported  by  the  striking  similarities  of  even  the 
most  complex  patterns  in  our  experiments  with  the  results  of 
the  simulations  in  open  systems  (for  example,  compare  our 
Fig.  7b  with  Fig.  12e  of  ref.  2),  which  represent  true  attractors. 

The  main  conclusions  are 

(1)  The  dynamical  phase  behavior  of  a  chain  of  BZ  drops  is  a 
function  of  coupling  strength.  As  coupling  strength  increases, 
the  following  sequence  is  observed  in  experiment  and  numer¬ 
ical  models  (Fig.  4c):  anti-phase,  bistable  anti-phase  and  in¬ 
phase/stationary  Turing  states,  in-phase. 

(2)  Malonic  acid  concentration  controls  the  coupling  between 
BZ  drops  in  oil  by  varying  the  balance  between  excitatory  and 
inhibitory  coupling.  Mechanistically,  malonic  acid  removes  the 
inhibitor;  therefore  decreasing  malonic  acid  increases  inhibitory 
coupling.  Drop  size  is  a  more  dominant  factor  in  coupling 
strength  than  the  length  of  the  oil  gap.  Theory  suggests  inhibitory 
coupling  strength  is  characterized  by  S  =  (PD)/(a(a  +  b)ke ff). 

(3)  Numerical  models  demonstrate  that  weak  coupling  solely 
through  the  inhibitors,  Br2  and  Br_,  produces  phase  repulsive 
coupling,  leading  to  anti-phase  synchrony,  while  coupling  solely  with 
activators,  HBr02  or  BrO/,  produces  attractive  phase  coupling, 
leading  to  in-phase  synchrony.  Strong  inhibitory  coupling  also 
produces  in-phase  synchrony;  essentially  the  drops  lose  their 
individual  identities  and  effectively  act  as  a  single  drop  in  the 
strong  coupling  limit.  The  transition  from  weak  to  strong 
coupling  is  marked  by  a  transition  from  out-of-phase  to  in-phase 
synchrony  and  occurs  at  S'  ~  1. 

(4)  Comparing  a  finite  element  model  to  experiment,  we 
conclude  that  Br2,  Br02#,  and  HBr02  exchange  between  drops 
with  partition  coefficients  PBr2  ^  2.5,  PBro2  ^  1,  and  PHBro2  *=  0.01; 
both  excitatory  and  inhibitory  coupling  need  to  be  included  in  the 
models  to  agree  with  the  experiments. 

(5)  The  point  oscillator  model  qualitatively  agrees  with 
experiment  and  a  realistic  3D  finite  element  model. 
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(6)  An  accurate  phase  model  was  constructed  for  the  case  of 
weak  coupling  and  anti-phase  attractors. 

(7)  Malonic  acid  concentration  decreases  with  time  in  the 
closed  system  of  microfluidic  drops.  In  contrast,  our  numerical 
models  consider  the  malonic  acid  concentration  to  be  constant. 

(8)  Malonic  acid  and  bromine  react  in  the  oil,  which  is  not 
accounted  for  in  our  numerical  models. 

(9)  Simple  numerical  models  account  for  the  majority  of  the 
observations. 
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